---
title: "COVID-19 Analysis"
author: "Kyle Fischer"
date: "13/12/2022"
output: html_document
---
Built with `r getRversion()`  

The purpose of this document is to analyse the full dataset on political ideology and reaction to the COVID19 pandemic.  

```{r setup, echo=F, warning=F, message=F}

library(tidyverse)
library(leaps)
library(MuMIn)
library(car)
library(Hmisc)
library(ordinal)
library(lavaan)
library(ggpubr)
library(lme4)
library(psych)
library(boot)
library(insight)
library(GPArotation)
library(gmodels)

# set ggplot theme
theme_set(theme_classic())
knitr::opts_chunk$set(echo = TRUE)

```

# Load Data & Remove Invalid Participants

```{r eval=F, echo=F, message=F, warning=F}

# Load data
  d <- read.csv("COVIDStudyData.csv")
  
# Perform Likert check calculations
  reverse <- function(x) 8 - x
  
  likertframe <- cbind(d$s2_sdo01, d$s2_sdo02, d$s2_sdo03, d$s2_sdo04, reverse(d$s2_sdo05r), reverse(d$s2_sdo06r), reverse(d$s2_sdo07r), 
                       reverse(d$s2_sdo08r), d$s2_sdo09, d$s2_sdo10, d$s2_sdo11, d$s2_sdo12, reverse(d$s2_sdo13r), reverse(d$s2_sdo14r), 
                       reverse(d$s2_sdo15r), reverse(d$s2_sdo16r), reverse(d$s2_rwa01r), d$s2_rwa02, reverse(d$s2_rwa03r), d$s2_rwa04, 
                       d$s2_rwa05, reverse(d$s2_rwa06r), reverse(d$s2_rwa07r), reverse(d$s2_rwa08r), reverse(d$s2_rwa09r), d$s2_rwa10, 
                       d$s2_rwa11, d$s2_rwa12, reverse(d$s2_rwa13r), d$s2_rwa14, d$s2_rwa15, reverse(d$s2_rwa16r), d$s2_rwa17, 
                       reverse(d$s2_rwa18r), reverse(d$s2_rwa19r), reverse(d$s2_rwa20r), d$s2_rwa21, d$s2_rwa22, d$s2_rwa23, 
                       reverse(d$s2_rwa24r), reverse(d$s2_rwa25r), d$s2_rwa26, reverse(d$s2_rwa27r), d$s2_rwa28, reverse(d$s2_rwa29r), 
                       d$s2_rwa30, d$s2_rwa31, d$s2_rwa32, reverse(d$s2_rwa33r), reverse(d$s2_rwa34r), d$s2_rwa35, reverse(d$s2_rwa36r), 
                       d$washingHands, d$sanitisingHands, d$waiveMedFees, d$grantPaidLeave, d$concernedForVulnerable, d$physDisOwn, 
                       d$facemask, d$banUSA, d$banChina, d$followGovRules, d$restrictMovement, d$restictEntryToUK, d$stronglyEnforceSD, 
                       d$severelyPunishViolators, reverse(d$forceStayHome_R), reverse(d$controlActions_R), d$bringInArmy, d$othersMustSD, 
                       d$dominantLeader, d$socMedShare, d$banPubGather, d$banRelGather, d$likelyavoidPubgather, d$likelyavoidRelgather,
                       d$spendMoneyOnVaccine, d$moreGovernmentResearch, d$impactEconomyWorseThanP, d$deathsEconomyWorseThanP,
                       d$government_response, d$likelyowndeath, d$likelyotherdeath, reverse(d$likelymeetotherR), d$unpleasantStayHome, 
                       d$airtravel, d$prestigeLeader, d$giveMoneyEconomy, d$stockUp)
  
   d <- d %>% mutate(s2_likert1 =  rowSums(likertframe == "1"), s2_likert2 =  rowSums(likertframe == "2"), 
                     s2_likert3 =  rowSums(likertframe == "3"), s2_likert4 =  rowSums(likertframe == "4"), 
                     s2_likert5 =  rowSums(likertframe == "5"), s2_likert6 =  rowSums(likertframe == "6"),
                     s2_likert7 =  rowSums(likertframe == "7"))

# Remove Invalid Participants
  d <- d %>%
    ## Study 1 exclusion criteria: 
        # TIME: less than 8 minutes/more than 90 minutes
        filter(S1_duration >= 480 & S1_duration <= 5400) %>%
        # SCALES: answer with the same response on more than 2/3rds of Likert scale questions
        filter(s1_badscales <= 40) %>%                                                                            
        # ATTENTION CHECKS: incorrectly answer at least 1 of 2 attention checks
        mutate(attentionFail1 = ifelse(S1_attention1 == "Strongly disagree", 0, 1),
               attentionFail2 = ifelse(S1_attention2 == "Slightly agree", 0, 1),
               numAttentionFails1 = attentionFail1 + attentionFail2) %>%
        filter(numAttentionFails1 %in% 0:1) %>%
        
    ## Study 2 exclusion criteria: 
        # TIME: less than 6 minutes/more than 60 minutes
        filter(S2_duration >= 360 & S2_duration <= 3600) %>%
        # SCALES: answer with the same response on more than 2/3rds of Likert scale questions
        filter(s2_likert1 <= 59 & s2_likert2 <= 59 & s2_likert3 <= 59 & s2_likert4 <= 59 & s2_likert5 <= 59 & 
               s2_likert6 <= 59 & s2_likert7 <= 59) %>%
        # ATTENTION CHECKS: incorrectly answer both of 2 attention checks
        mutate(attentionFail3 = ifelse(S2_attention1 == "Slightly agree", 0, 1),
               attentionFail4 = ifelse(S2_attention2 == "Slightly agree", 0, 1),
               numAttentionFails2 = attentionFail3 + attentionFail4) %>%
        filter(numAttentionFails2 == 0)

# Cleanup
  rm(reverse, likertframe)
  
```

# Create Additional Variables

```{r eval=F, echo=F, message=F, warning=F}

# Means                                                                              
  d <- d %>% # Ideology means 
             mutate(s1_meanSDO             = rowMeans(select(., starts_with("s1_sdo"))),
                    s1_meanRWA             = rowMeans(select(., starts_with("s1_rwa"))),
                    s2_meanSDO             = rowMeans(select(., starts_with("s2_sdo"))),
                    s2_meanRWA             = rowMeans(select(., starts_with("s2_rwa"))),
                    s3_meanEmpathConcern   = rowMeans(cbind(d$EC01, d$EC03r, d$EC05, d$EC07r, d$EC09r, d$EC10, d$EC12))) %>%
             # COVID19
             mutate(s2_meanOtherRegard     = rowMeans(cbind(d$waiveMedFees, d$grantPaidLeave, d$concernedForVulnerable)),
                    s2_meanAuth            = rowMeans(cbind(d$followGovRules, d$othersMustSD, d$restrictMovement, d$controlActions_R, 
                                                            d$forceStayHome_R, d$stronglyEnforceSD, d$banChina, d$banUSA, 
                                                            d$restictEntryToUK, d$severelyPunishViolators, d$bringInArmy, d$moAnger, 
                                                            d$moDisgust, d$moContempt, d$moOutrage, d$moKTreated)))

# Exposure variable
  d <- d %>% mutate(exposureAll = ifelse(diagnosed == "Yes" | 
                                         knowDiagnosedPerson == "Yes" | 
                                         suspectHad == "Yes" |
                                         medTreatDelayedSelf == "Yes" | 
                                         medTreatDelayedOther == "Yes" | 
                                         financialImpact == "Yes", 1, 0))
  
# Covariates
  d <- d %>% mutate(ethnicityCat = ifelse(ethnicityS == "White", 0, 1),
                    polCat       = ifelse(voting %in% c("Labour Party", "Green Party"), "LEFT", ifelse(voting %in% c("Conservatives",
                                   "UK Independence Party"), "RIGHT", ifelse(pol == "Left", "LEFT", ifelse(pol == "Right", "RIGHT", 
                                   ifelse(pol == "Centre", "CENTRE", NA))))))
# Remove polCat for participants whose political responses don't match
  d <- d %>% mutate(polCat = ifelse((pol == "Left" & voting == "Conservatives") |
                                    (pol == "Right" & (voting == "Green Party" | voting == "Labour Party")),
                                    NA, polCat))
# Set polCat indicator variables  
  d <- d %>% mutate(polCatR      = ifelse(is.na(polCat) == 1, NA, ifelse(polCat == "RIGHT", 1, 0)),
                    polCatC      = ifelse(is.na(polCat) == 1, NA, ifelse(polCat == "CENTRE", 1, 0)))

# Change in SDO/RWA
  d <- d %>% mutate(changeSDO = s2_meanSDO - s1_meanSDO, changeRWA = s2_meanRWA - s1_meanRWA)

```

# Data Statistics

```{r echo=F, message=F, warning=F}

# Main Text
# All 427 participants in both wave 1 & 2
  describe(d$age)$min
  describe(d$age)$max
  describe(d$age)$mean
  table(d$sex, useNA = "ifany") ## Female is 0
  table(d$ethnicityCat, useNA = "ifany") ## "White" is 0
  table(d$polCat, useNA = "ifany")

# Participants from wave 3 (remove s3 participants With completion time outside of 3 deviations around the mean)
  s3minDur <- ifelse(mean(d$S3_duration, na.rm = TRUE) - (3 * sd(d$S3_duration, na.rm = TRUE)) < 0, 0, 
                   mean(d$S3_duration, na.rm = TRUE) - (3 * sd(d$S3_duration, na.rm = TRUE)))
  s3maxDur <- mean(d$S3_duration, na.rm = TRUE) + (3 * sd(d$S3_duration, na.rm = TRUE))
  dEC      <- d %>% filter(S3_duration > s3minDur & S3_duration < s3maxDur)
  
  describe(dEC$age)$min
  describe(dEC$age)$max
  describe(dEC$age)$mean
  table(dEC$sex, useNA = "ifany")   ## Female is 0
  table(dEC$ethnicityCat, useNA = "ifany") ## "White" is 0
  table(dEC$polCat, useNA = "ifany")
  
# Missing responses
  # Wave 2 - no political affiliation
  dPol <- d %>% filter(is.na(polCat) == 0) #380
  # Wave 2 - no political affiliation & no sliders (worry/MO questions)
  dPolNS <- dPol %>% filter(is.na(worryPandemic) == 0) #374
  # Wave 3 no political affiliation
  dECPol <- dEC %>% filter(is.na(polCat) == 0) #312
  # Wave 3 no political affiliation & no sliders
  dECPolNS <- dECPol %>% filter(is.na(worryPandemic) == 0) #308

# Supplementary Materials Section 1
  # Dropouts wave 1 to 2: See additional RMD file
  
  # Dropouts wave 2 to 3: 
  d <- d %>% mutate(wave3 = ifelse(is.na(S3_duration), 0, ifelse(S3_duration > s3minDur & S3_duration < s3maxDur, 1, 0)))
  print("SDO")
  t.test(s1_meanSDO ~ wave3, data = select(d, "wave3", "s1_meanSDO"), var.equal = TRUE)   #if non-significant then the means are equal
  var.test(s1_meanSDO ~ wave3, data = select(d, "wave3", "s1_meanSDO"))                  #if non-significant then variances equal
  print("RWA")
  t.test(s1_meanRWA ~ wave3, data = select(d, "wave3", "s1_meanRWA"), var.equal = TRUE)  #if non-significant then the means are equal
  var.test(s1_meanRWA ~ wave3, data = select(d, "wave3", "s1_meanRWA"))                 #if non-significant then variances equal
  print("political affiliation")  
  CrossTable(d$wave3, d$polCat, fisher = TRUE, chisq = FALSE, expected = TRUE,
  sresid = TRUE, format = "SPSS", prop.c = FALSE, prop.t = FALSE)
  print("ses")
  t.test(ses ~ wave3, data = select(d, "wave3", "ses"), var.equal = TRUE)  #if non-significant then the means are equal
  var.test(ses ~ wave3, data = select(d, "wave3", "ses"))                 #if non-significant then variances equal
  print("sex")  
  CrossTable(d$wave3, d$sex, fisher = TRUE, chisq = FALSE, expected = TRUE,
  sresid = TRUE, format = "SPSS", prop.c = FALSE, prop.t = FALSE)
  print("age")
  t.test(age ~ wave3, data = select(d, "wave3", "age"), var.equal = TRUE)  #if non-significant then the means are equal
  var.test(age ~ wave3, data = select(d, "wave3", "age"))                  #if non-significant then variances equal
  print("race")  
  CrossTable(d$wave3, d$ethnicityCat, fisher = TRUE, chisq = FALSE, expected = TRUE,
  sresid = TRUE, format = "SPSS", prop.c = FALSE, prop.t = FALSE)

```

# Reliability: SDO, RWA, Empathic concern

```{r eval=F, echo=F, message=F, warning=F}

# Calculate reliability
  print("SDO S1")
  alpha(cov(data.matrix(d[,c("s1_sdo01", "s1_sdo02", "s1_sdo03", "s1_sdo04", "s1_sdo05r", "s1_sdo06r", "s1_sdo07r", "s1_sdo08r", "s1_sdo09", 
                             "s1_sdo10", "s1_sdo11", "s1_sdo12", "s1_sdo13r", "s1_sdo14r", "s1_sdo15r", "s1_sdo16r")])))
  
  print("RWA S1")
  alpha(cov(data.matrix(d[,c("s1_rwa01r", "s1_rwa02", "s1_rwa03r", "s1_rwa04", "s1_rwa05", "s1_rwa06r", "s1_rwa07r", "s1_rwa08r", "s1_rwa09r", 
                             "s1_rwa10", "s1_rwa11", "s1_rwa12", "s1_rwa13r", "s1_rwa14", "s1_rwa15", "s1_rwa16r", "s1_rwa17", "s1_rwa18r", 
                             "s1_rwa19r", "s1_rwa20r", "s1_rwa21", "s1_rwa22", "s1_rwa23", "s1_rwa24r", "s1_rwa25r", "s1_rwa26", "s1_rwa27r", 
                             "s1_rwa28", "s1_rwa29r", "s1_rwa30", "s1_rwa31", "s1_rwa32", "s1_rwa33r", "s1_rwa34r", "s1_rwa35", "s1_rwa36r")])))
  
  print("empathic concern")
  alpha(cov(data.matrix(dEC[,c("EC01", "EC03r", "EC05", "EC07r", "EC09r", "EC10", "EC12")])))

```

# Fit measures for Structural Equation Modelling - Confirmatory Factor Analyses

```{r eval=F, echo=F, message=F, warning=F}

# RMSEA (Root Mean Square Error of Approximation): A parsimony-adjusted index. Values closer to 0 represent a good fit. Cut-off: < 0.08
# SRMR (Standardized Root Mean Square Residual):   The square-root of the difference between the residuals of the sample covariance matrix 
#                                                  and the hypothesized model. Cut-off: < 0.08

# SDO
  print("Pre-pandemic SDO")
  model    <- 'SDO1 =~ s1_sdo01 + s1_sdo02 + s1_sdo03 + s1_sdo04 + s1_sdo05r + s1_sdo06r + s1_sdo07r + s1_sdo08r + 
                       s1_sdo09 + s1_sdo10 + s1_sdo11 + s1_sdo12 + s1_sdo13r + s1_sdo14r + s1_sdo15r + s1_sdo16r'
  cfa_SDO1 <- cfa(model, data = d, ordered = c("s1_sdo01", "s1_sdo02", "s1_sdo03", "s1_sdo04", "s1_sdo05r", "s1_sdo06r", "s1_sdo07r", "s1_sdo08r", 
                                               "s1_sdo09", "s1_sdo10", "s1_sdo11", "s1_sdo12", "s1_sdo13r", "s1_sdo14r", "s1_sdo15r", "s1_sdo16r"))
  fitMeasures(cfa_SDO1)["rmsea"] %>% round(2)
  fitMeasures(cfa_SDO1)["srmr"] %>% round(2)

# RWA
  print("Pre-pandemic RWA")
  model <- 'RWA1 =~ s1_rwa01r + s1_rwa02 + s1_rwa03r + s1_rwa04 + s1_rwa05 + s1_rwa06r + s1_rwa07r + s1_rwa08r + s1_rwa09r + s1_rwa10 + 
                    s1_rwa11 + s1_rwa12 + s1_rwa13r + s1_rwa14 + s1_rwa15 + s1_rwa16r + s1_rwa17 + s1_rwa18r + s1_rwa19r + s1_rwa20r + 
                    s1_rwa21 + s1_rwa22 + s1_rwa23 + s1_rwa24r + s1_rwa25r + s1_rwa26 + s1_rwa27r + s1_rwa28 + s1_rwa29r + s1_rwa30 + 
                    s1_rwa31 + s1_rwa32 + s1_rwa33r + s1_rwa34r + s1_rwa35 + s1_rwa36r'
  cfa_RWA1 <- cfa(model, data = d, ordered = c("s1_rwa01r", "s1_rwa02", "s1_rwa03r", "s1_rwa04", "s1_rwa05", "s1_rwa06r", "s1_rwa07r", "s1_rwa08r", 
                                               "s1_rwa09r", "s1_rwa10", "s1_rwa11", "s1_rwa12", "s1_rwa13r", "s1_rwa14", "s1_rwa15", "s1_rwa16r", 
                                               "s1_rwa17", "s1_rwa18r", "s1_rwa19r", "s1_rwa20r", "s1_rwa21", "s1_rwa22", "s1_rwa23", "s1_rwa24r", 
                                               "s1_rwa25r", "s1_rwa26", "s1_rwa27r", "s1_rwa28", "s1_rwa29r", "s1_rwa30", "s1_rwa31", "s1_rwa32", 
                                               "s1_rwa33r", "s1_rwa34r", "s1_rwa35", "s1_rwa36r"))
  fitMeasures(cfa_RWA1)["rmsea"] %>% round(2)
  fitMeasures(cfa_RWA1)["srmr"] %>% round(2)

## Can SEM performance be improved by removing scale items?

# SDO short version
  print("Pre-pandemic SDO short version")
  model    <- 'SDO1 =~ s1_sdo03 + s1_sdo04 + s1_sdo05r + s1_sdo06r + s1_sdo11 + s1_sdo12 + s1_sdo13r + s1_sdo14r'
  cfa_SDO1 <- cfa(model, data = d, ordered = c("s1_sdo03", "s1_sdo04", "s1_sdo05r", "s1_sdo06r", "s1_sdo11", "s1_sdo12", "s1_sdo13r", "s1_sdo14r"))
  fitMeasures(cfa_SDO1)["rmsea"] %>% round(2)
  fitMeasures(cfa_SDO1)["srmr"] %>% round(2)

# RWA short version
  print("Pre-pandemic RWA short version")
  model <- 'RWA1 =~ s1_rwa01r + s1_rwa02 + s1_rwa03r + s1_rwa04 + s1_rwa05 + s1_rwa06r + s1_rwa13r + s1_rwa14 + s1_rwa15 + s1_rwa16r + s1_rwa17 + 
            s1_rwa18r + s1_rwa25r + s1_rwa26 + s1_rwa27r + s1_rwa28 + s1_rwa29r + s1_rwa30'
  cfa_RWA1 <- cfa(model, data = d, ordered = c("s1_rwa01r", "s1_rwa02", "s1_rwa03r", "s1_rwa04", "s1_rwa05", "s1_rwa06r", "s1_rwa13r", "s1_rwa14", 
                                               "s1_rwa15", "s1_rwa16r", "s1_rwa17", "s1_rwa18r", "s1_rwa25r", "s1_rwa26", "s1_rwa27r", "s1_rwa28", 
                                               "s1_rwa29r", "s1_rwa30"))
  fitMeasures(cfa_RWA1)["rmsea"] %>% round(2)
  fitMeasures(cfa_RWA1)["srmr"] %>% round(2)
  rm(cfa_RWA1, cfa_SDO1) 
  
# Model fit was not adequate for RWA, and many items would need to be removed to adequately improve fit.
  
```

# Principal Component Analyses

```{r eval=F, echo=F, message=F, warning=F}

# Create dataset:
  PCAData <- d %>% filter(is.na(worryPandemic) == 0) %>% 
                   select (followGovRules,           #“It is important to follow the UK government's rules regarding COVID-19”
                           othersMustSD,             #“Because of COVID-19, it is very important that others take physical distancing very seriously and 
                                                     #limit all social contact”
                           restrictMovement,         #“I support government measures to restrict the movement of UK citizens to limit the spread of COVID-19”
                           controlActions_R,         #“It makes me angry that the government would tell me where I can go and what I can do, even when 
                                                     #there is a crisis..."
                           forceStayHome_R,          #“I am upset at the thought that my government would force people to stay home against their will”
                           stronglyEnforceSD,        #“It is vital right now that the government strongly enforces social distancing measures”
                           banChina,                 #“All citizens of China should be banned from entering the UK while the COVID-19 pandemic is ongoing”
                           banUSA,                   #“All citizens of the USA should be banned from entering the UK while the COVID-19 pandemic is ongoing”     
                           restictEntryToUK,         #“Strict entry restrictions should be imposed at all borders while the COVID-19 pandemic is ongoing”
                           severelyPunishViolators,  #“I want my government to severely punish those who violate orders to stay home”
                           bringInArmy,              #“The army should be mobilised to enforce quarantines and rules regarding COVID-19”
                           moAnger,                  #"To what extent does K’s behaviour make you feel anger?"
                           moDisgust,                #"To what extent does K’s behaviour make you feel disgust?"
                           moContempt,               #"To what extent does K’s behaviour make you feel contempt?"
                           moOutrage,                #"To what extent does K’s behaviour make you feel outrage?"
                           moKTreated)               #“Finally, how do you think K should be treated?”

# Check correlations: 

  # 1. Check correlation matrix:
  PCAMatrix <- cor(PCAData)
  # Low correlations: Variables that only have a small number of correlations greater than .3
  view(colSums(PCAMatrix > 0.3))
  # Multicollinearity: Any correlations greater than .9
  view(colSums(PCAMatrix > 0.9)) 
  
  # 2. Check Kaiser–Meyer–Olkin measure: 
  # just >.5 : barely acceptable, .5-.7 : mediocre, .7-.8 : good, .8-.9 : great, >.9 : superb (Hutcheson & Sofroniou, 1999)
  # If 0 factor analysis likely inappropriate, 1 factor analysis should yield distinct and reliable factors
  # Should be above the bare minimum of .5 for all variables as well as overall (preferably higher)
  KMO(PCAData)
  
  # 3. Run Bartlett’s test of sphericity on the correlation matrix: 
  # If significant, factor analysis is appropriate.
  cortest.bartlett(PCAData)          
  
  # 4. Check determinant: 
  # If greater than 0.00001, determinant is not problematic.
  det(PCAMatrix) 
  
# Determine number of components to use:
  
  # 1. Use a scree plot to find point of inflection:
  PCA1 <- principal(PCAData, nfactors = 16, rotate = "none")
  PCA1
  plot(PCA1$values, type = "b")
  fa.parallel(PCAData, fm = 'minres', fa = 'fa')       

  # 2. Rerun PCA, with the number of components, and:
  # a. Check Kaiser’s criterion - that:
  #    - If there are fewer than 30 variables, all communalities (h2) are greater than .7, or 
  #    - If the sample size exceeds 250, the average communality (h2) is greater than .6.
  # b. Check all SS loadings are greater than 0.7 (SS loadings = eigenvalues, Proportion Var = as a proportion of the no of factors (i.e. 3) - gives 
  # the variance explained)
  # c. Check fit based upon off diagonal values: > 0.95 considered indicator of good fit
  PCA2 <- principal(PCAData, nfactors = 4, rotate = "none")
  PCA2
  plot(PCA2$values, type="b")

  # 3. Check residuals:
  # Most must be less than 0.05, if proportion more than 50% have grounds for concern
  PCARes <- factor.residuals(PCAMatrix, PCA2$loadings)
  PCARes <- as.matrix(PCARes[upper.tri(PCARes)])
  PCALargeRes <- abs(PCARes) > 0.05
  sum(PCALargeRes) / nrow(PCARes)
  # If root-mean-square residual high (> 0.08) consider extracting more factors
  sqrt(mean(PCARes^2))
  # Distribution should be approximately normal
  hist(PCARes)
  
# Print Factor Loading (Pattern) Matrix  (theory suggests factors might correlate so use oblimin)
  # Display only loadings above .3 and sort items by the size of loadings
  PCAFLM <- principal(PCAData, nfactors = 4, rotate = "oblimin")
  print.psych(PCAFLM, cut = 0.3, sort = TRUE)
  
  print.psych(PCAFLM, sort = TRUE)
  
# Check reliability: 
  # Raw α: overall & individually: items greater than "overall" may need to be deleted from the scale to improve reliability
  # rdrop: if less than about .3, that item does not correlate very well with the scale overall
  # Frequencies (final table): percentage of each response to each items - make sure that everyone is not giving the same response
  #  Usually if everyone (or almost everyone) gives the same response on an item it will have poor reliability stats.
  alpha(PCAData[,c(1,2,3,4,5,6)])     # Support for lockdown rules
  alpha(PCAData[,c(12,13,14,15,16)])  # Moral emotions towards rule-breakers
  alpha(PCAData[,c(7,8,9)])           # Support for strict border control
  alpha(PCAData[,c(10,11)])           # Support for severe government enforcement of rules

# Add PC data to dataset
  d <- d %>% mutate(PC1SupportLockdownRules       = rowMeans(cbind(d$followGovRules, d$othersMustSD, d$restrictMovement, d$controlActions_R, 
                                                                   d$forceStayHome_R, d$stronglyEnforceSD)),
                    PC2MoralEmotions              = rowMeans(cbind(d$moAnger, d$moDisgust, d$moContempt, d$moOutrage, d$moKTreated)),
                    PC3SupportStrictBorderControl = rowMeans(cbind(d$banChina, d$banUSA, d$restictEntryToUK)),
                    PC4SupportSevereEnforcement   = rowMeans(cbind(d$severelyPunishViolators, d$bringInArmy)))
  
# Save final data file
  write.csv(file = "FinalDataFile.csv", d, row.names = FALSE)
   
rm(PCAData, PCAMatrix, PCA1, PCA2, PCARes, PCALargeRes, PCAFLM, PCASLM)
  
```

# Set up datasets

```{r eval=F, echo=F, message=F, warning=F}

# Standardized variables for GLMs
  dZ <- d %>% select(id, s1_meanSDO, s1_meanRWA, s2_meanSDO, s2_meanRWA, s3_meanEmpathConcern, s2_meanOtherRegard, s2_meanAuth,  
                     PC1SupportLockdownRules, PC2MoralEmotions, PC3SupportStrictBorderControl, PC4SupportSevereEnforcement, changeRWA, changeSDO, 
                     worryPandemic, worrySelfGetSick, worryFamilyGetSick, worryOthersGetSick, worryFinancial, worryEconomy, physDisOwn, facemask, 
                     socMedShare, age, sex, ses, ethnicityCat, pol, polCat, polCatR, polCatC, exposureAll, S3_duration, voting)  
  
  dZ <- dZ %>% mutate(s1_meanSDO                   = (s1_meanSDO - mean(dZ$s1_meanSDO, na.rm = TRUE)) / sd(dZ$s1_meanSDO, na.rm = TRUE),
                      s1_meanRWA                   = (s1_meanRWA - mean(dZ$s1_meanRWA, na.rm = TRUE)) / sd(dZ$s1_meanRWA, na.rm = TRUE),
                      s2_meanSDO                   = (s2_meanSDO - mean(dZ$s2_meanSDO, na.rm = TRUE)) / sd(dZ$s2_meanSDO, na.rm = TRUE),
                      s2_meanRWA                   = (s2_meanRWA - mean(dZ$s2_meanRWA, na.rm = TRUE)) / sd(dZ$s2_meanRWA, na.rm = TRUE),
                      s3_meanEmpathConcern         = (s3_meanEmpathConcern - mean(dZ$s3_meanEmpathConcern, na.rm = TRUE)) / sd(dZ$s3_meanEmpathConcern, 
                                                                                                                               na.rm = TRUE),
                      s2_meanOtherRegard            = (s2_meanOtherRegard - mean(dZ$s2_meanOtherRegard, na.rm = TRUE)) / sd(dZ$s2_meanOtherRegard, na.rm = TRUE),
                      s2_meanAuth                   = (s2_meanAuth - mean(dZ$s2_meanAuth, na.rm = TRUE)) / sd(dZ$s2_meanAuth, na.rm = TRUE),
                      PC1SupportLockdownRules       = (PC1SupportLockdownRules - mean(dZ$PC1SupportLockdownRules, na.rm = TRUE)) / 
                                                       sd(dZ$PC1SupportLockdownRules, na.rm = TRUE),
                      PC2MoralEmotions              = (PC2MoralEmotions - mean(dZ$PC2MoralEmotions, na.rm = TRUE)) / sd(dZ$PC2MoralEmotions, na.rm = TRUE),
                      PC3SupportStrictBorderControl = (PC3SupportStrictBorderControl - mean(dZ$PC3SupportStrictBorderControl, na.rm = TRUE)) / 
                                                       sd(dZ$PC3SupportStrictBorderControl, na.rm = TRUE),
                      PC4SupportSevereEnforcement   = (PC4SupportSevereEnforcement - mean(dZ$PC4SupportSevereEnforcement, na.rm = TRUE)) / 
                                                       sd(dZ$PC4SupportSevereEnforcement, na.rm = TRUE),
                      changeRWA                     = (changeRWA - mean(dZ$changeRWA, na.rm = TRUE)) / sd(dZ$changeRWA, na.rm = TRUE),
                      changeSDO                     = (changeSDO - mean(dZ$changeSDO, na.rm = TRUE)) / sd(dZ$changeSDO, na.rm = TRUE),
                      worryPandemic                 = (worryPandemic - mean(dZ$worryPandemic, na.rm = TRUE)) / sd(dZ$worryPandemic, na.rm = TRUE),
                      worrySelfGetSick              = (worrySelfGetSick - mean(dZ$worrySelfGetSick,  na.rm = TRUE)) / sd(dZ$worrySelfGetSick, na.rm = TRUE),
                      worryFamilyGetSick            = (worryFamilyGetSick - mean(dZ$worryFamilyGetSick, na.rm = TRUE)) / sd(dZ$worryFamilyGetSick, na.rm = TRUE),
                      worryOthersGetSick            = (worryOthersGetSick - mean(dZ$worryOthersGetSick, na.rm = TRUE)) / sd(dZ$worryOthersGetSick, na.rm = TRUE),
                      worryFinancial                = (worryFinancial - mean(dZ$worryFinancial, na.rm = TRUE)) / sd(dZ$worryFinancial, na.rm = TRUE),
                      worryEconomy                  = (worryEconomy - mean(dZ$worryEconomy, na.rm = TRUE)) / sd(dZ$worryEconomy, na.rm = TRUE),
                      physDisOwn                    = (physDisOwn - mean(dZ$physDisOwn, na.rm = TRUE)) / sd(dZ$physDisOwn, na.rm = TRUE),
                      facemask                      = (facemask - mean(dZ$facemask, na.rm = TRUE)) / sd(dZ$facemask, na.rm = TRUE),
                      socMedShare                   = (socMedShare - mean(dZ$socMedShare, na.rm = TRUE)) / sd(dZ$socMedShare, na.rm = TRUE),
                      age                           = (age - mean(dZ$age, na.rm = TRUE)) / sd(dZ$age, na.rm = TRUE),
                      ses                           = (ses - mean(dZ$ses, na.rm = TRUE)) / sd(dZ$ses, na.rm = TRUE))

# Empathic concern (remove s3 participants With completion time outside of 3 deviations around the mean)
  dZEC   <- dZ %>% filter(S3_duration > s3minDur & S3_duration < s3maxDur)
  rm(s3minDur, s3maxDur)
      
```

# Descriptive stats for all relevant variables

```{r}

# Supplementary Materials - Section 3

# Means & standard deviations
  rbind(cbind(Variable = 'Socioeconomic status', Mean = round(mean(d$ses),2), SD = round(sd(d$ses),2)),
        cbind('Empathic concern', round(mean(dEC$s3_meanEmpathConcern),2), round(sd(dEC$s3_meanEmpathConcern),2)),
        cbind('SDO (before pandemic)', round(mean(d$s1_meanSDO),2), round(sd(d$s1_meanSDO),2)),
        cbind('SDO (during pandemic)', round(mean(d$s1_meanRWA),2), round(sd(d$s1_meanRWA),2)),
        cbind('RWA (before pandemic)', round(mean(d$s2_meanSDO),2), round(sd(d$s2_meanSDO),2)),
        cbind('RWA (during pandemic)', round(mean(d$s2_meanRWA),2), round(sd(d$s2_meanRWA),2)),
        cbind('Cooperative/other-regarding responses to COVID-19', round(mean(d$s2_meanOtherRegard),2), round(sd(d$s2_meanOtherRegard),2)),
        cbind('Conformist/norm-enforcing responses to COVID-19', round(mean(d$s2_meanAuth, na.rm = TRUE),2), round(sd(d$s2_meanAuth, na.rm = TRUE),2)),
        cbind('Support for lockdown rules (PC1)', round(mean(d$PC1SupportLockdownRules),2), round(sd(d$PC1SupportLockdownRules),2)),
        cbind('Moral emotions towards rule breaker (PC2)', round(mean(d$PC2MoralEmotions, na.rm = TRUE),2), round(sd(d$PC2MoralEmotions, na.rm = TRUE),2)),
        cbind('Support strict border control (PC3)', round(mean(d$PC3SupportStrictBorderControl),2), round(sd(d$PC3SupportStrictBorderControl),2)),
        cbind('Support severe enforcement (PC4)', round(mean(d$PC4SupportSevereEnforcement, na.rm = TRUE),2), round(sd(d$PC4SupportSevereEnforcement),2)), 
        cbind('Concern about COVID-19', round(mean(d$worryPandemic, na.rm = TRUE),2), round(sd(d$worryPandemic, na.rm = TRUE),2)),
        cbind('Concern about own health', round(mean(d$worrySelfGetSick, na.rm = TRUE),2), round(sd(d$worrySelfGetSick, na.rm = TRUE),2)),
        cbind('Concern about familiar others’ health', round(mean(d$worryFamilyGetSick, na.rm = TRUE),2), round(sd(d$worryFamilyGetSick, na.rm = TRUE),2)),
        cbind('Concern about unfamiliar others’ health', round(mean(d$worryOthersGetSick, na.rm = TRUE),2), round(sd(d$worryOthersGetSick, na.rm = TRUE),2)),
        cbind('Concern about own finances', round(mean(d$worryFinancial, na.rm = TRUE),2), round(sd(d$worryFinancial, na.rm = TRUE),2)),  
        cbind('Concern about economy', round(mean(d$worryEconomy, na.rm = TRUE),2), round(sd(d$worryEconomy, na.rm = TRUE),2)))    

# Histograms  
  histogram(d$ses)
  histogram(dEC$s3_meanEmpathConcern)
  histogram(d$s1_meanSDO)
  histogram(d$s1_meanRWA)
  histogram(d$s2_meanOtherRegard)
  histogram(d$s2_meanAuth)
  histogram(d$PC1SupportLockdownRules)
  histogram(d$PC2MoralEmotions)
  histogram(d$PC3SupportStrictBorderControl)
  histogram(d$PC4SupportSevereEnforcement)
  histogram(d$worryPandemic)
  histogram(d$worrySelfGetSick)
  histogram(d$worryFamilyGetSick)
  histogram(d$worryOthersGetSick)
  histogram(d$worryFinancial)
  histogram(d$worryEconomy)

```

# Correlation Matrix

```{r eval=F, echo=F, message=F, warning=F}

# Correlation matrix
  cor_data <- d %>% select(s1_meanSDO, s1_meanRWA,
                           s3_meanEmpathConcern, 
                           s2_meanOtherRegard, s2_meanAuth, 
                           PC1SupportLockdownRules, PC2MoralEmotions, PC3SupportStrictBorderControl, PC4SupportSevereEnforcement, 
                           worryPandemic, worrySelfGetSick, worryFamilyGetSick, worryOthersGetSick, worryFinancial, worryEconomy,
                           age, ses)

  cor_matrix <- rcorr(as.matrix(cor_data))
  cor_matrix_x <- rbind(cor_matrix$r, cor_matrix$P)

  write.csv(cor_matrix_x, file = "Correlation matrix.csv")

  rm(cor_data, cor_matrix, cor_matrix_x)

```

# Set up regression functions

```{r eval=F, echo=F, message=F, warning=F}

# General linear models (GLM) functions
  modelGLM <- function(m)
    {
        s <- summary(m)
        data.frame(Model = as.character(s$call)[2], 
                   n     = n_obs(m),
                   Fstat = round(s$fstatistic[1],2), 
                   DF1   = round(s$fstatistic[2],2), 
                   DF2   = round(s$fstatistic[3],2), 
                   p     = round(pf(s$fstatistic[1], s$fstatistic[2], s$fstatistic[3], lower.tail = FALSE),3), 
                   Rsq   = round(s$r.squared,2))
    }

  coefficientsGLM <- function(m, name)
    {
        s <- summary(m)
        data.frame('Name' = c(rep(name, times = attributes(logLik(m))$df - 1)),
                   c(data.frame(IV = rownames(as.data.frame(round(s$coefficients,6))))),
                   c(data.frame(round(s$coefficients,6))),
                   round(confint(m),2))
    } 
  
  bootReg <- function (formula, data, indices)
    {
        d <- data [indices,]
        fit <- lm(formula, data = d)
        return(coef(fit))
    }

  bootCI <- function(model, formula, data)
    {
       bootResults <- boot(statistic = bootReg, formula = formula, data = data, R = 2000)
       bootCoef    <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
       s <- summary(model)                          
       for(i in 1:(attributes(logLik(model))$df) - 1)
       {
          a <- boot.ci(bootResults, type = "bca", index = i)  
          bootCoef <- rbind(bootCoef, cbind(name = rownames(as.data.frame(s$coefficients))[i], BCa.lower = as.numeric(a$bca[4]), BCa.upper = as.numeric(a$bca[5])))    
       }
       return(bootCoef)
    }

  plotGLM <- function(data, name, width, height, font) 
    {
                    glmPlot <- ggplot(data,                
                                      aes(x     = forcats::fct_rev(factor(Name)),
                                          y     = Estimate,
                                          color = IV,
                                          shape = IV)) +  
                      
                               geom_point(size  = 4,
                                          position = pd) +
                          
                               geom_errorbar(aes(ymin  = lowerCI,
                                                 ymax  = upperCI),
                                                 width = 0.1,
                                                 size  = 0.5,
                                                 position = pd) +
                               theme_bw() +
                               theme(axis.text.y  = element_text(size = font),
                                     axis.title.x = element_text(size = font),
                                     legend.text  = element_text(size = font)) +
                      
                               geom_hline(yintercept = 0, color = "black", size = 1) +        
                               ylab("Standardised Beta Estimates with Bootstrapped Confidence Intervals") + xlab("") + 
                               coord_flip()  +
                               scale_shape_manual(values = c(16, 18)) + 
                               scale_colour_manual(values = c("red", "blue"))

                       ggsave(name, plot = glmPlot, device = "jpeg", width = width, height = height, dpi = 300) 
                       glmPlot
    }

  plotGLMA <- function(data, name, width, height, font) 
    {
                    glmPlot <- ggplot(data,                
                                      aes(x     = forcats::fct_rev(factor(Name)),
                                          y     = Estimate,
                                          color = IV,
                                          shape = IV)) +  
                      
                               geom_point(size  = 4,
                                          position = pd) +
                          
                               geom_errorbar(aes(ymin  = lowerCI,
                                                 ymax  = upperCI),
                                                 width = 0.1,
                                                 size  = 0.5,
                                                 position = pd) +
                               theme_bw() +
                               theme(axis.text.y  = element_text(size = font),
                                     axis.title.x = element_text(size = font),
                                     legend.text  = element_text(size = font)) +
                      
                               geom_hline(yintercept = 0, color = "black", size = 1) +        
                               ylab("Standardised Beta Estimates with Bootstrapped Confidence Intervals") + xlab("") + 
                               coord_flip()  +
                               scale_shape_manual(values = c(16, 15, 18, 17)) + #RWAC, RWA + 17 is the triangle, 18 is the diamond
                               scale_colour_manual(values = c("red", "red", "blue", "blue"))

                       ggsave(name, plot = glmPlot, device = "jpeg", width = width, height = height, dpi = 300) 
                       glmPlot
    }

  plotGLMB <- function(data, name, width, height, font) 
    {
                    glmPlot <- ggplot(data,                
                                      aes(x     = forcats::fct_rev(factor(Name)),
                                          y     = Estimate,
                                          color = IV,
                                          shape = IV)) +  
                      
                               geom_point(size  = 4,
                                          position = pd) +
                          
                               geom_errorbar(aes(ymin  = lowerCI,
                                                 ymax  = upperCI),
                                                 width = 0.1,
                                                 size  = 0.5,
                                                 position = pd) +
                               theme_bw() +
                               theme(axis.text.y  = element_text(size = font),
                                     axis.title.x = element_text(size = font),
                                     legend.text  = element_text(size = font)) +
                      
                               geom_hline(yintercept = 0, color = "black", size = 1) +        
                               ylab("Standardised Beta Estimates with Bootstrapped Confidence Intervals") + xlab("") + 
                               coord_flip()  +
                               scale_shape_manual(values = c(8, 4, 15, 16, 17, 18)) + 
                               scale_colour_manual(values = c("green", "green", "red", "red", "blue", "blue"))

                       ggsave(name, plot = glmPlot, device = "jpeg", width = width, height = height, dpi = 300) 
                       glmPlot
    }

 plotGLMC <- function(data, name, width, height, font) 
    {
                    glmPlot <- ggplot(data,                
                                      aes(x     = forcats::fct_rev(factor(Name)),
                                          y     = Estimate,
                                          color = IV,
                                          shape = IV)) +  
                      
                               geom_point(size  = 4,
                                          position = pd) +
                          
                               geom_errorbar(aes(ymin  = lowerCI,
                                                 ymax  = upperCI),
                                                 width = 0.1,
                                                 size  = 0.5,
                                                 position = pd) +
                               coord_cartesian(xlim = c(-1, 1)) +
                               theme_bw() +
                               theme(axis.text.y  = element_text(size = font),
                                     axis.title.x = element_text(size = font),
                                     legend.text  = element_text(size = font)) +
                      
                               geom_hline(yintercept = 0, color = "black", size = 1) +        
                               ylab("Standardised Beta Estimates with Bootstrapped Confidence Intervals") + xlab("") + 
                               coord_flip()  +
                               scale_shape_manual(values = c(16, 15, 18, 17)) + #RWAC, RWA + 17 is the triangle, 18 is the diamond
                               scale_colour_manual(values = c("red", "red", "blue", "blue"))

                       ggsave(name, plot = glmPlot, device = "jpeg", width = width, height = height, dpi = 300) 
                       glmPlot
    }

     
  pd = position_dodge(0.7)  #How much to jitter the points on plots
  pd2 = position_dodge(1.0)  #How much to jitter the points on plots

```

# Prediction 1: Cooperative and Conformist Responses
## 1a. SDO and RWA and mean of cooperative/mean of conformist items

```{r eval=F, echo=F, message=F, warning=F}

# Setup
  names    <- vector(mode = "list", length = 4)
  models   <- vector(mode = "list", length = 4)
  formulas <- vector(mode = "list", length = 4)
  P1.models       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  P1.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  P1.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())

# Models
  names[[1]]    <- "1. Cooperative/other-regarding"
  models[[1]]   <- lm(s2_meanOtherRegard ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[1]] <- s2_meanOtherRegard ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[2]]    <- "2. Conformist/norm-enforcing" 
  models[[2]]   <- lm(s2_meanAuth ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[2]] <- s2_meanAuth ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
 
# Results
  for(i in 1:2) 
  {
      P1.models         <- rbind(P1.models, modelGLM(models[[i]]))
      P1.coefficients   <- rbind(P1.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      P1.bootcoef       <- rbind(P1.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  P1.coefficients <- cbind(P1.coefficients, data.frame(BCa.lower = as.numeric(as.character(P1.bootcoef$BCa.lower)), 
                                                       BCa.upper = as.numeric(as.character(P1.bootcoef$BCa.upper))))

  # View results
  P1.results <- P1.coefficients %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  P1.results <- P1.results %>% mutate(IV = ifelse(IV == "s1_meanRWA", "RWA", "SDO")) %>%
                               mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  P1.results <- P1.results %>% select(Name, Estimate, lowerCI, upperCI, IV)

# Plot
  plotGLM (P1.results, "glmPlotCooperationConformity.png", 10, 3, 10)

##########################################################################################################################################################
  
# Supplementary materials: models without political affiliation:

# Temp tables
  temp.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  temp.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
 
# Models
  names[[3]]    <- "1. Cooperative/other-regarding (excl. political affiliation)"
  models[[3]]   <- lm(s2_meanOtherRegard ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[3]] <- s2_meanOtherRegard ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[4]]    <- "2. Conformist/norm-enforcing (excl. political affiliation)" 
  models[[4]]   <- lm(s2_meanAuth ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[4]] <- s2_meanAuth ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll

# Results
  for(i in 3:4) 
  {
      P1.models           <- rbind(P1.models, modelGLM(models[[i]]))
      temp.coefficients   <- rbind(temp.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      temp.bootcoef       <- rbind(temp.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  temp.coefficients <- cbind(temp.coefficients, data.frame(BCa.lower = as.numeric(as.character(temp.bootcoef$BCa.lower)), 
                                                       BCa.upper = as.numeric(as.character(temp.bootcoef$BCa.upper))))

  # View results
  temp.results <- temp.coefficients %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  temp.results <- temp.results %>% 
                    mutate(IV = ifelse(IV == "s1_meanRWA", "RWA excl for political affiliation", "SDO excl for political affiliation")) %>%
                    mutate(Name = ifelse(Name == "1. Cooperative/other-regarding (excl. political affiliation)", 
                                        "1. Cooperative/other-regarding", "2. Conformist/norm-enforcing")) %>%
                    mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  temp.results <- temp.results %>% select(Name, Estimate, lowerCI, upperCI, IV)
  
  P1.coefficients.full <- rbind(P1.coefficients, temp.coefficients)
  P1.results.full <- rbind(P1.results, temp.results)
  
# Plot
  plotGLMA (P1.results.full, "glmPlotCooperationConformity (supp).png", 10, 3, 10)

  # view(P1.models %>% mutate_if(is.numeric, round, digits = 4))
  # view(P1.coefficients.full %>% mutate_if(is.numeric, round, digits = 4))
  
rm(plotData, names, models, formulas, i, P1.bootcoef, temp.models, temp.coefficients, temp.bootcoef)

```

## 1b. SDO and RWA and components of conformist items

```{r eval=F, echo=F, message=F, warning=F}

# Setup
  names    <- vector(mode = "list", length = 8)
  models   <- vector(mode = "list", length = 8)
  formulas <- vector(mode = "list", length = 8)
  P1.modelsb       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  P1.coefficientsb <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  P1.bootcoefb     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())

# Models
  names[[1]]    <- "1. Support lockdown rules (PC1)"
  models[[1]]   <- lm(PC1SupportLockdownRules ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[1]] <- PC1SupportLockdownRules ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[2]]    <- "2. Moral emotions towards rule breaker (PC2)"
  models[[2]]   <- lm(PC2MoralEmotions ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[2]] <- PC2MoralEmotions ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[3]]    <- "3. Support strict border control (PC3)"
  models[[3]]   <- lm(PC3SupportStrictBorderControl ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[3]] <- PC3SupportStrictBorderControl ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[4]]    <- "4. Support severe enforcement (PC4)"
  models[[4]]   <- lm(PC4SupportSevereEnforcement ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[4]] <- PC4SupportSevereEnforcement ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
    
# Results
  for(i in 1:4) 
  {
      P1.modelsb         <- rbind(P1.modelsb, modelGLM(models[[i]]))
      P1.coefficientsb   <- rbind(P1.coefficientsb, coefficientsGLM(models[[i]], names[[i]]))
      P1.bootcoefb       <- rbind(P1.bootcoefb, bootCI(models[[i]], formulas[[i]], dZ))
  }
  P1.coefficientsb <- cbind(P1.coefficientsb, data.frame(BCa.lower = as.numeric(as.character(P1.bootcoefb$BCa.lower)), 
                                                       BCa.upper = as.numeric(as.character(P1.bootcoefb$BCa.upper))))
  
# View results
  P1.resultsb <- P1.coefficientsb %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  P1.resultsb <- P1.resultsb %>% mutate(IV = ifelse(IV == "s1_meanRWA", "RWA", "SDO")) %>%
                                 mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  P1.resultsb <- P1.resultsb %>% select(Name, Estimate, lowerCI, upperCI, IV)

  # Plot
  plotGLM (P1.resultsb, "glmPlotCooperationConformityB.png", 10, 4, 10)

##########################################################################################################################################################
  
# Supplementary materials: models without political affiliation:

# Temp tables
  temp.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  temp.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
  
# Models
  names[[5]]    <- "1. Support lockdown rules (PC1) (excl. political affiliation)"
  models[[5]]   <- lm(PC1SupportLockdownRules ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[5]] <- PC1SupportLockdownRules ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[6]]    <- "2. Moral emotions towards rule breaker (PC2) (excl. political affiliation)"
  models[[6]]   <- lm(PC2MoralEmotions ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[6]] <- PC2MoralEmotions ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[7]]    <- "3. Support strict border control (PC3) (excl. political affiliation)"
  models[[7]]   <- lm(PC3SupportStrictBorderControl ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[7]] <- PC3SupportStrictBorderControl ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[8]]    <- "4. Support severe enforcement (PC4) (excl. political affiliation)"
  models[[8]]   <- lm(PC4SupportSevereEnforcement ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[8]]  <- PC4SupportSevereEnforcement ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll

# Results
  for(i in 5:8) 
  {
      P1.modelsb          <- rbind(P1.modelsb, modelGLM(models[[i]]))
      temp.coefficients   <- rbind(temp.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      temp.bootcoef       <- rbind(temp.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  temp.coefficients <- cbind(temp.coefficients, data.frame(BCa.lower = as.numeric(as.character(temp.bootcoef$BCa.lower)), 
                                                       BCa.upper = as.numeric(as.character(temp.bootcoef$BCa.upper))))

# View results
  temp.results <- temp.coefficients %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  temp.results <- temp.results %>% 
                    mutate(IV = ifelse(IV == "s1_meanRWA", "RWA excl for political affiliation", "SDO excl for political affiliation")) %>%
                    mutate(Name = ifelse(Name == "1. Support lockdown rules (PC1) (excl. political affiliation)", "1. Support lockdown rules (PC1)",
                                   ifelse(Name == "2. Moral emotions towards rule breaker (PC2) (excl. political affiliation)", 
                                                        "2. Moral emotions towards rule breaker (PC2)",
                                    ifelse(Name == "3. Support strict border control (PC3) (excl. political affiliation)",
                                                        "3. Support strict border control (PC3)", "4. Support severe enforcement (PC4)")))) %>%
                    mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  temp.results <- temp.results %>% select(Name, Estimate, lowerCI, upperCI, IV)
  
  P1.coefficientsb.full <- rbind(P1.coefficientsb, temp.coefficients)
  P1.resultsb.full <- rbind(P1.resultsb, temp.results)

# Plot  
  plotGLMA (P1.resultsb.full, "glmPlotCooperationConformityB (supp).png", 10, 4, 10)

  # view(P1.modelsb %>% mutate_if(is.numeric, round, digits = 4))
  # view(P1.coefficientsb.full %>% mutate_if(is.numeric, round, digits = 4))

rm(plotData, names, models, formulas, i, P1.bootcoef, temp.models, temp.coefficients, temp.bootcoef)
   
```

## 1c. Empathic concern

```{r eval=F, echo=F, message=F, warning=F}

# Setup
  names    <- vector(mode = "list", length = 8)
  models   <- vector(mode = "list", length = 8)
  formulas <- vector(mode = "list", length = 8)
  P1.modelsc       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  P1.coefficientsc <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  P1.bootcoefc     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())                                                     

# Models  
  ## 3a: RWA & EC with political affiliation
  names[[1]]    <- "EC: Support lockdown rules (PC1)"
  models[[1]]   <- lm(PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZEC)
  formulas[[1]] <- PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[2]]    <- "EC: Moral emotions towards rule breaker (PC2)"
  models[[2]]   <- lm(PC2MoralEmotions ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZEC)
  formulas[[2]] <- PC2MoralEmotions ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[3]]    <- "EC: Support strict border control (PC3)"
  models[[3]]   <- lm(PC3SupportStrictBorderControl ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZEC)
  formulas[[3]] <- PC3SupportStrictBorderControl ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[4]]    <- "EC: Support severe enforcement (PC4)"
  models[[4]]   <- lm(PC4SupportSevereEnforcement ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZEC)
  formulas[[4]] <- PC4SupportSevereEnforcement ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  ## 3b: RWA & EC without political affiliation
  names[[5]]    <- "EC: Support lockdown rules (PC1) (excl. political affiliation)"
  models[[5]]   <- lm(PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
  formulas[[5]] <- PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[6]]    <- "EC: Moral emotions towards rule breaker (PC2) (excl. political affiliation)"
  models[[6]]   <- lm(PC2MoralEmotions ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
  formulas[[6]] <- PC2MoralEmotions ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[7]]    <- "EC: Support strict border control (PC3) (excl. political affiliation)"
  models[[7]]   <- lm(PC3SupportStrictBorderControl ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
  formulas[[7]] <- PC3SupportStrictBorderControl ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[8]]    <- "EC: Support severe enforcement (PC4) (excl. political affiliation)"
  models[[8]]   <- lm(PC4SupportSevereEnforcement ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
  formulas[[8]] <- PC4SupportSevereEnforcement ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll  
  
# Results
  for(i in 1:8) 
  {
      P1.modelsc         <- rbind(P1.modelsc, modelGLM(models[[i]]))
      P1.coefficientsc   <- rbind(P1.coefficientsc, coefficientsGLM(models[[i]], names[[i]]))
      P1.bootcoefc       <- rbind(P1.bootcoefc, bootCI(models[[i]], formulas[[i]], dZ))
  }
  P1.coefficientsc.full <- cbind(P1.coefficientsc, data.frame(BCa.lower = as.numeric(as.character(P1.bootcoefc$BCa.lower)), 
                                                         BCa.upper = as.numeric(as.character(P1.bootcoefc$BCa.upper))))
  
  # view(P1.modelsc %>% mutate_if(is.numeric, round, digits = 4))
  # view(P1.coefficientsc.full %>% mutate_if(is.numeric, round, digits = 4))

rm(names, models, formulas, i, P1.bootcoefc)

```

## 1d. SDO/empathic concern mediation analyses

```{r eval=F, echo=F, message=F, warning=F}

#Is the relationship between empathy and PC1 mediated by SDO?

library(mediation)
#fitM <- lm(M ~ X, data = d)
#fitY <- lm(Y ~ X + M, data = d)

summary(lm(s1_meanSDO ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + polCatR + polCatC + exposureAll, data = dZEC))

########## empathy -> SDO -> support for lockdown rules

model_mediator <- lm(s1_meanSDO ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + polCatR + polCatC + exposureAll, data = dZEC)
model_outcome  <- lm(PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + polCatR + polCatC + exposureAll, data = dZEC)

?mediate
mediation_result <- mediate(model_mediator, model_outcome, boot = TRUE, sims = 1000, treat = "s3_meanEmpathConcern", mediator = "s1_meanSDO", dropobs = TRUE)

summary(mediation_result)
plot(mediation_result)

# no political affiliation

summary(lm(s1_meanSDO ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC))

model_mediator <- lm(s1_meanSDO ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
model_outcome  <- lm(PC1SupportLockdownRules ~ s3_meanEmpathConcern + s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)

?mediate
mediation_result <- mediate(model_mediator, model_outcome, boot = TRUE, sims = 1000, treat = "s3_meanEmpathConcern", mediator = "s1_meanSDO", dropobs = TRUE)

summary(mediation_result)
plot(mediation_result)

########## empathy -> SDO -> moral emotions
# no political affiliation

model_mediator <- lm(s1_meanSDO ~ s3_meanEmpathConcern + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)
model_outcome  <- lm(PC2MoralEmotions ~ s3_meanEmpathConcern + s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZEC)

?mediate
mediation_result <- mediate(model_mediator, model_outcome, boot = TRUE, sims = 1000, treat = "s3_meanEmpathConcern", mediator = "s1_meanSDO", dropobs = TRUE)

summary(mediation_result)
plot(mediation_result)

```

# Prediction 2: Threats

```{r eval=F, echo=F, message=F, warning=F}

# Set up
  names    <- vector(mode = "list", length = 12)
  models   <- vector(mode = "list", length = 12)
  formulas <- vector(mode = "list", length = 12)
  P3.models       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  P3.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  P3.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
  
# Models
  names[[1]]   <- "1. Worry about pandemic"
  models[[1]]  <- lm(worryPandemic ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[1]] <- worryPandemic ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[2]]   <- "2. Worry about getting sick"
  models[[2]]  <- lm(worrySelfGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[2]] <- worrySelfGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[3]]   <- "3. Worry about familiar others getting sick"
  models[[3]]  <- lm(worryFamilyGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[3]] <- worryFamilyGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[4]]   <- "4. Worry about unfamiliar others getting sick"
  models[[4]]  <- lm(worryOthersGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[4]] <- worryOthersGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[5]]   <- "5. Worry about personal finances"
  models[[5]]  <- lm(worryFinancial ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[5]] <- worryFinancial ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[6]]   <- "6. Worry about economy"
  models[[6]]  <- lm(worryEconomy ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[6]] <- worryEconomy ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC

# Results
  for(i in 1:6) 
  {
      P3.models       <- rbind(P3.models, modelGLM(models[[i]]))
      P3.coefficients <- rbind(P3.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      P3.bootcoef     <- rbind(P3.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  P3.coefficients <- cbind(P3.coefficients, data.frame(BCa.lower = as.numeric(as.character(P3.bootcoef$BCa.lower)), 
                                                       BCa.upper = as.numeric(as.character(P3.bootcoef$BCa.upper)))) 
  
  P3.results <- P3.coefficients %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  P3.results <- P3.results %>% mutate(IV = ifelse(IV == "s1_meanRWA", "RWA", "SDO")) %>% 
                               mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  P3.results <- P3.results %>% dplyr::select(Name, Estimate, lowerCI, upperCI, IV)

# Plot
  plotGLM (P3.results, "glmPlotThreat.png", 10, 6, 10)

##########################################################################################################################################################
  
# Supplementary materials: models without political affiliation:

# Temp tables
  temp.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  temp.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
  
# Models
  names[[7]]   <- "1. Worry about pandemic (excl. political affiliation)"
  models[[7]]  <- lm(worryPandemic ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[7]] <- worryPandemic ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[8]]   <- "2. Worry about getting sick (excl. political affiliation)"
  models[[8]]  <- lm(worrySelfGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[8]] <- worrySelfGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[9]]   <- "3. Worry about familiar others getting sick (excl. political affiliation)"
  models[[9]]  <- lm(worryFamilyGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[9]] <- worryFamilyGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[10]]   <- "4. Worry about unfamiliar others getting sick (excl. political affiliation)"
  models[[10]]  <- lm(worryOthersGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[10]] <- worryOthersGetSick ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[11]]   <- "5. Worry about personal finances (excl. political affiliation)"
  models[[11]]  <- lm(worryFinancial ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[11]] <- worryFinancial ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll
  names[[12]]   <- "6. Worry about economy (excl. political affiliation)"
  models[[12]]  <- lm(worryEconomy ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[12]] <- worryEconomy ~ s1_meanSDO + s1_meanRWA + age + sex + ethnicityCat + ses + exposureAll

# Results
  for(i in 7:12) 
  {
      P3.models         <- rbind(P3.models, modelGLM(models[[i]]))
      temp.coefficients <- rbind(temp.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      temp.bootcoef     <- rbind(temp.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  temp.coefficients <- cbind(temp.coefficients, data.frame(BCa.lower = as.numeric(as.character(temp.bootcoef$BCa.lower)), 
                                                           BCa.upper = as.numeric(as.character(temp.bootcoef$BCa.upper)))) 
  
  temp.results <- temp.coefficients %>% filter(IV == "s1_meanRWA" | IV == "s1_meanSDO")
  temp.results <- temp.results %>% mutate(IV = ifelse(IV == "s1_meanRWA", "RWA excl for political affiliation", "SDO excl for political affiliation"))
  temp.results <- temp.results %>% mutate(Name = ifelse(Name == "1. Worry about pandemic (excl. political affiliation)" , "1. Worry about pandemic", 
                                                  ifelse(Name == "2. Worry about getting sick (excl. political affiliation)" , "2. Worry about getting sick",
                                                    ifelse(Name == "3. Worry about familiar others getting sick (excl. political affiliation)", 
                                                                    "3. Worry about familiar others getting sick",
                                                      ifelse(Name == "4. Worry about unfamiliar others getting sick (excl. political affiliation)", 
                                                                      "4. Worry about unfamiliar others getting sick",
                                                        ifelse(Name == "5. Worry about personal finances (excl. political affiliation)", "5. Worry about personal finances",
                                                                        "6. Worry about economy"))))))
  temp.results <- temp.results %>% mutate(lowerCI = BCa.lower, upperCI = BCa.upper)
  temp.results <- temp.results %>% dplyr::select(Name, Estimate, lowerCI, upperCI, IV)

  P3.coefficients.full <- rbind(P3.coefficients, temp.coefficients)
  P3.results.full <- rbind(P3.results, temp.results)
    
# Plot
  plotGLMC (P3.results.full, "glmPlotThreat (supp).png", 10, 6, 10)
      
# View results
  # view(P3.models %>% mutate_if(is.numeric, round, digits = 4))
  # view(P3.coefficients.full %>% mutate_if(is.numeric, round, digits = 4))
 
rm(names, models, modelsZ, formulas, datasets, i, P3.bootcoef, temp.bootcoef, temp.results)

```

# Prediction 3: Increase in response to COVID
## 3a. Was there an increase in RWA/SDO?

```{r eval=F, echo=F, message=F, warning=F}

# RWA
  ## Before
  mean(d$s1_meanRWA, na.rm = TRUE); sqrt(var(d$s1_meanRWA, na.rm = TRUE)/length(d$s1_meanRWA))
  ## During  
  mean(d$s2_meanRWA, na.rm = TRUE); sqrt(var(d$s2_meanRWA, na.rm = TRUE)/length(d$s2_meanRWA))
   
  ## t-test
  t <- t.test(d$s2_meanRWA, d$s1_meanRWA, paired = TRUE, alternative = "greater")
  t
  round((t$statistic[[1]]) ^ 2 / ((t$statistic[[1]]) ^ 2 + (t$parameter[[1]])), 3) # r-squared
     
# SDO   
  ## SDO Before
  mean(d$s1_meanSDO, na.rm = TRUE); sqrt(var(d$s1_meanSDO, na.rm = TRUE)/length(d$s1_meanSDO))
  ## SDO During
  mean(d$s2_meanSDO, na.rm = TRUE); sqrt(var(d$s2_meanSDO, na.rm = TRUE)/length(d$s2_meanSDO))
      
  ## t-test
  t <- t.test(d$s2_meanSDO, d$s1_meanSDO, paired = TRUE, alternative = "two.sided") 
  t
  ## r-squared:
  round((t$statistic[[1]]) ^ 2 / ((t$statistic[[1]]) ^ 2 + (t$parameter[[1]])), 3)
  
``` 
   
## 3b. Was the increase in RWA driven by perceived and real pandemic-related threats?

```{r eval=F, echo=F, message=F, warning=F}

# Set up
  names    <- vector(mode = "list", length = 4)
  models   <- vector(mode = "list", length = 4)
  formulas <- vector(mode = "list", length = 4)
  P2.models       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  P2.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  P2.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
    
# Models
  names[[1]]   <- "Worrying about the pandemic"
  models[[1]]  <- lm(changeRWA ~ worryPandemic + age + sex + ethnicityCat + ses + polCatR + polCatC, data = dZ)
  formulas[[1]] <- changeRWA ~ worryPandemic + age + sex + ethnicityCat + ses + polCatR + polCatC
  names[[2]]   <- "Exposure to COVID-19"
  models[[2]]  <- lm(changeRWA ~ exposureAll + age + sex + ethnicityCat + ses + polCatR + polCatC, data = dZ)
  formulas[[2]] <- changeRWA ~ exposureAll + age + sex + ethnicityCat + ses + polCatR + polCatC
  
  names[[3]]   <- "Worrying about the pandemic (excl. political affiliation)"
  models[[3]]  <- lm(changeRWA ~ worryPandemic + age + sex + ethnicityCat + ses, data = dZ)
  formulas[[3]] <- changeRWA ~ worryPandemic + age + sex + ethnicityCat + ses
  names[[4]]   <- "Exposure to COVID-19 (excl. political affiliation)"
  models[[4]]  <- lm(changeRWA ~ exposureAll + age + sex + ethnicityCat + ses, data = dZ)
  formulas[[4]] <- changeRWA ~ exposureAll + age + sex + ethnicityCat + ses

# Results
  for(i in 1:4) 
  {
      P2.models       <- rbind(P2.models, modelGLM(models[[i]]))
      P2.coefficients <- rbind(P2.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      P2.bootcoef     <- rbind(P2.bootcoef, bootCI(models[[i]], formulas[[i]], dZEC))
  }
  P2.coefficients.full <- cbind(P2.coefficients, data.frame(BCa.lower = as.numeric(as.character(P2.bootcoef$BCa.lower)), 
                                                            BCa.upper = as.numeric(as.character(P2.bootcoef$BCa.upper))))

# View results
  #view(P2.models %>% mutate_if(is.numeric, round, digits = 4))
  #view(P2.coefficients %>% mutate_if(is.numeric, round, digits = 4))
  
rm(names, models, formulas, i, P2.bootcoef)

```

# Interaction Analyses
## With political affiliation

```{r eval=F, echo=F, message=F, warning=F}

# Set up
  names    <- vector(mode = "list", length = 12)
  models   <- vector(mode = "list", length = 12)
  formulas <- vector(mode = "list", length = 12)
  IA.models       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  IA.coefficients <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  IA.bootcoef     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
  
# Models
  ## Means
  names[[1]]    <- "Interaction: Mean cooperative/other-regarding"
  models[[1]]   <- lm(s2_meanOtherRegard ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[1]] <- s2_meanOtherRegard ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[2]]    <- "Interaction: Mean conformist/norm-enforcing"
  models[[2]]   <- lm(s2_meanAuth ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[2]] <- s2_meanAuth ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  ## PCs
  names[[3]]    <- "Interaction: Support lockdown rules (PC1)"
  models[[3]]   <- lm(PC1SupportLockdownRules ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[3]] <- PC1SupportLockdownRules ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[4]]    <- "Interaction: Moral emotions towards rule breaker (PC2)"
  models[[4]]   <- lm(PC2MoralEmotions ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[4]] <- PC2MoralEmotions ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[5]]    <- "Interaction: Support strict border control (PC3)"
  models[[5]]   <- lm(PC3SupportStrictBorderControl ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[5]] <- PC3SupportStrictBorderControl ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[6]]    <- "Interaction: Support severe enforcement (PC4)"
  models[[6]]   <- lm(PC4SupportSevereEnforcement ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[6]] <- PC4SupportSevereEnforcement ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  ## Concern
  names[[7]]    <- "Interaction: Worrying about pandemic"
  models[[7]]   <- lm(worryPandemic ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[7]] <- worryPandemic ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[8]]    <- "Interaction: Worry about getting sick"
  models[[8]]   <- lm(worrySelfGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[8]] <- worrySelfGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[9]]    <- "Interaction: Worry about familiar others getting sick"
  models[[9]]   <- lm(worryFamilyGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[9]] <- worryFamilyGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[10]]    <- "Interaction: Worry about unfamiliar others getting sick"
  models[[10]]   <- lm(worryOthersGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[10]] <- worryOthersGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[11]]    <- "Interaction: Worry about personal finances"
  models[[11]]   <- lm(worryFinancial ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[11]] <- worryFinancial ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  names[[12]]    <- "Interaction: Worry about economy"
  models[[12]]   <- lm(worryEconomy ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC, data = dZ)
  formulas[[12]] <- worryEconomy ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll + polCatR + polCatC
  
# Results
  for(i in 1:12) 
  {
      IA.models       <- rbind(IA.models, modelGLM(models[[i]]))
      IA.coefficients <- rbind(IA.coefficients, coefficientsGLM(models[[i]], names[[i]]))
      IA.bootcoef     <- rbind(IA.bootcoef, bootCI(models[[i]], formulas[[i]], dZ))
  }
  IA.coefficients.full <- cbind(IA.coefficients, data.frame(BCa.lower = as.numeric(as.character(IA.bootcoef$BCa.lower)), 
                                                            BCa.upper = as.numeric(as.character(IA.bootcoef$BCa.upper)))) 
  
# View results
  #view(IA.models %>% mutate_if(is.numeric, round, digits = 4))
  #view(IA.coefficients.full %>% mutate_if(is.numeric, round, digits = 4))
  
  rm(names, models, modelsZ, formulas, datasets, i, IA.bootcoef)

```

## Without political affiliation

```{r eval=F, echo=F, message=F, warning=F}

# Set up
  names    <- vector(mode = "list", length = 12)
  models   <- vector(mode = "list", length = 12)
  formulas <- vector(mode = "list", length = 12)
  IA.models.noPA       <- data.frame(Model=character(), n=numeric(), Fstat=numeric(), DF1=numeric(), DF2=numeric(), p=numeric(),Rsq=numeric(), Rsqa=numeric())
  IA.coefficients.noPA <- data.frame(DV=character(), IV=character(), estimate=numeric(), se=numeric(), tValue=numeric(), pvalue=numeric(),
                              'est2.5'=numeric(), 'est97.5'=numeric(), estimateSTD=numeric(), seSTD=numeric(), tValueSTD=numeric(), pvalueSTD=numeric(),
                              'est2.5STD'=numeric(), 'est97.5STD'=numeric(), stringsAsFactors=FALSE)
  IA.bootcoef.noPA     <- data.frame(name=character(), BCa.lower = numeric(), BCa.upper = numeric())
  
# Models
  ## Means
  names[[1]]    <- "Interaction: Mean cooperative/other-regarding (excl. political affiliation)"
  models[[1]]   <- lm(s2_meanOtherRegard ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[1]] <- s2_meanOtherRegard ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[2]]    <- "Interaction: Mean conformist/norm-enforcing (excl. political affiliation)"
  models[[2]]   <- lm(s2_meanAuth ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[2]] <- s2_meanAuth ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  ## PCs  
  names[[3]]    <- "Interaction: Support lockdown rules (PC1) (excl. political affiliation)"
  models[[3]]   <- lm(PC1SupportLockdownRules ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[3]] <- PC1SupportLockdownRules ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[4]]    <- "Interaction: Moral emotions towards rule breaker (PC2) (excl. political affiliation)"
  models[[4]]   <- lm(PC2MoralEmotions ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[4]] <- PC2MoralEmotions ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[5]]    <- "Interaction: Support strict border control (PC3) (excl. political affiliation)"
  models[[5]]   <- lm(PC3SupportStrictBorderControl ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[5]] <- PC3SupportStrictBorderControl ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[6]]    <- "Interaction: Support severe enforcement (PC4) (excl. political affiliation)"
  models[[6]]   <- lm(PC4SupportSevereEnforcement ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[6]] <- PC4SupportSevereEnforcement ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  ## Concern
  names[[7]]    <- "Interaction: Worrying about pandemic (excl. political affiliation)"
  models[[7]]   <- lm(worryPandemic ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[7]] <- worryPandemic ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[8]]    <- "Interaction: Worry about getting sick (excl. political affiliation)"
  models[[8]]   <- lm(worrySelfGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[8]] <- worrySelfGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[9]]    <- "Interaction: Worry about familiar others getting sick (excl. political affiliation)"
  models[[9]]   <- lm(worryFamilyGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[9]] <- worryFamilyGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[10]]    <- "Interaction: Worry about unfamiliar others getting sick (excl. political affiliation)"
  models[[10]]   <- lm(worryOthersGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[10]] <- worryOthersGetSick ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[11]]    <- "Interaction: Worry about personal finances (excl. political affiliation)"
  models[[11]]   <- lm(worryFinancial ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[11]] <- worryFinancial ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  names[[12]]    <- "Interaction: Worry about economy (excl. political affiliation)"
  models[[12]]   <- lm(worryEconomy ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll, data = dZ)
  formulas[[12]] <- worryEconomy ~ s1_meanRWA*s1_meanSDO + age + sex + ethnicityCat + ses + exposureAll
  
# Results
  for(i in 1:12) 
  {
      IA.models.noPA       <- rbind(IA.models.noPA, modelGLM(models[[i]]))
      IA.coefficients.noPA <- rbind(IA.coefficients.noPA, coefficientsGLM(models[[i]], names[[i]]))
      IA.bootcoef.noPA     <- rbind(IA.bootcoef.noPA, bootCI(models[[i]], formulas[[i]], dZ))
  }
  IA.coefficients.noPA.full <- cbind(IA.coefficients.noPA, data.frame(BCa.lower = as.numeric(as.character(IA.bootcoef.noPA$BCa.lower)), 
                                                                      BCa.upper = as.numeric(as.character(IA.bootcoef.noPA$BCa.upper)))) 
  
# View results
  # view(IA.models.noPA %>% mutate_if(is.numeric, round, digits = 4))
  # view(IA.coefficients.noPA.full %>% mutate_if(is.numeric, round, digits = 4))
  
rm(names, models, modelsZ, formulas, datasets, i, IA.bootcoef.noPA)

```

### Save GLM Results

```{r eval=F, echo=F, message=F, warning=F}

# Section 4 Models
GLM.models <- rbind(P1.models, P1.modelsb, P1.modelsc, P3.models, P2.models, IA.models, IA.models.noPA)
GLM.coefficients <- rbind(P1.coefficients.full, P1.coefficientsb.full, P1.coefficientsc.full, P3.coefficients.full, P2.coefficients.full, IA.coefficients.full, IA.coefficients.noPA.full)

write.csv(file = "GLM - Model Stats.csv", GLM.models, row.names = FALSE)
write.csv(file = "GLM - Full Results.csv", GLM.coefficients, row.names = FALSE)

```
